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Abstract 

Multitime correlation functions provide useful probes for the ensembles of trajectories underlying 
the stochastic dynamics of complex systems. These can be obtained by measuring their optical 
response to sequences of ultrashort optical pulse. Using the continuous time random walk model 
for spectral diffusion, we analyze the signatures of anomalous relaxation in two-dimensional four 
wave mixing signals. Different models which share the same two point joint probability distribution 
show markedly different lineshapes and may be distinguished. Aging random walks corresponding 
to waiting time distributions with diverging first moment show dependence of 2D lineshapes on 
initial observation time, which persist for long times. 
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I. INTRODUCTION 



Simple relaxation theories break down when the relaxation is non-exponential and as- 
sumes, for example, a stretched-exponential or an algebraic form. Such anomalous re- 
laxation has been observed in numerous physical systems ranging from sirigle molecules 
and quantum do^ spectral diffusion in fluorescence blinking trajectories 
protein fol ding M, |8j, charge- carrier transport, geophysical processes, and in economics 



a, 



10 
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12 . Il3| . Stochastic dynamics can be fully described by ensembles of trajec- 



tories of collective variables [ij]. Statistical analysis of stochastic trajectories results in a 
hierarchy of multipoint correlation functions which carry increasing levels of information. 
Two-point correlation functions provide the simplest measure of fluctuations and the most 
common evidence for anomalous relaxation. They are the easiest to sample experimentally 
and to predict theoretically. However they do not uniquely characterize the system. Many 
models can be constructed that have the same two point correlations but very different 
higher order correlation functions. Anomalous dynamics implies that many timescales are 
relevant. These may represent various dynamical variables or metastable configurations in 
polymers or glassy systems 15|, ll6|, 113]. Treating all relevant variables explicitly is not always 



possible. Some calculations only include directly accessible variables (such as the transition 
frequency in spectral diffusion) [18] and use a master equation for their probability densities; 
all other variables are projected out and represented through memory functions. The long 
time memories characteristic of anomalous relaxation are not compatible with the ordinary 
Markovian a ppr oximation which assumes fast memory loss. The master equations derived 



in this case 
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20| are thus limited to two point correlation functions and do not 
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carry 



221. 



enough information to describe the multipoint correlation and response functions 

Several practical strategies may be employed towards the simulation of multipoint corre- 
lation functions. One option is to use Markovian master equations with a large number of 
collective variables. Another possibility is to assume harmonic (Gaussian) processes which 



are exactly solvable 23|, . All information is then contained in the spectral density, which 
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25H. A different class of solvable mod- 



27|, which assume the erasure of all 



may be tailored to give long-tailed correlations 
els are continuous time random walks (CTRW) 
memory (renewal) when the relevant dynamical variables are changed (jumps). They por- 
tray the dynamics as a generalized random walk with a distributed waiting time or length 
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for stochastic jumps between various states. Memory enters this model solely through the 
time t elapsed from last renewal time. Anomalous behavior is observed when the waiting 
time distribution function (WTDF) ip{t) for the next jump has long tails. We have recently 
proposed that lineshapes in coherent multidimensional optical spectroscopy may be used to 
probe anomalous multipoint correlation functions [28]. Algebraic singularities at transition 
frequencies and power-law cross-peak dynamics were predicted in the two-dimensional op- 
tical response of a two-level chromophore to three laser pulses whose frequency undergoes 
a stochastic two state jump continuous time random walk with a power-law waiting time 
density function iplt) ~ t~"~^. In this paper we present more detailed simulations for this 
model and further demonstrate how it may be used to probe aging effects in systems that 
never equilibrate. Frequency- domain signals such as linear absorption are ill defined in ag- 
ing systems since they depend on the measurement time window. 2DCS is a time-domain 
technique that uses ultrashort pulses. Such signals should provide unambiguous signatures 
for aging, since all delay times are fully controlled. Different models may be distinguished 
by higher order nonlinear techniques. 

We shall focus on two classes of WTDF which lead to anomalous spectral lineshapes. We 
assume asymptotic algebraic decay ijjit) ~ which shows significant deviations from 



normal relaxation for < a < 2 2^ 
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32]. For 1 < a < 2 stationary ensembles may 



be described by a proper choice of initial condition, which implies a special WTDF for the 
first jump tp'{t) which represents how the system was prepared. The anomalous multipoint 
correlations observed in fluorescence traces of conformation dynamics of flavin proteins 25 ] 
showed symmetries due to microscopic reversibility typical for stationary processes. 

For < a < 1, stationary ensembles can not be constructed. System properties neces- 
sarily depend on the time elapsed from the initial preparation even when it is very long. 
This phenomena is known as aging. Such random walks show fractal behavior related to 
^^evy stable distributions, which generalize the Gaussian distributions of ordinary diffusion 
291 ]. This case is fundamentally more complicated than 1 < a < 2: such random walks 
are nonergodic [33], their time and ensemble averages may differ 3^], and special sample 
preparation for each run of the experiment is needed. Signatures of aging were observed in 
fluorescence blinking of single CdSe quantum dots, with a ~ 0.5 Isl. I35I . This is in agree- 



ment with the Sparre- Andersen theorem 
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38| which states that the first passage time 



of random walk with any symmetric distribution of jump lengths (including Levy flights) 
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has a universal asymptotic ~ t~^/^ decay. The origin of these long tailed WTDF is not fully 
understood. 

Environment dynamics affects spectral lineshapes through modulations of the transition 
frequencies. However, extracting the fluctuation timescales from absorption lineshapes is 
not always possible and may require additional assumptions and the introduction of speciflc 
models. Nonlinear spectroscopies can distinguish between nonequivalent dynamical models 
whose linear response is identical. In two-Dimensional Correlation Spectroscopy (2008)0, 
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41| the system is subjected to three femtosecond laser pulses (Fig 1). The flrst pulse 



creates a coherence between the ground state and an excited state. Time evolution (free 
induction decay) during the flrst interval ti is related to the absorption lineshape by a Fourier 
transform. The second pulse erases the coherence, bringing the molecule to the ground or an 
excited state population. The transition frequency continues to change by interaction with 
the environment during the second interval t2- Finally a coherence is again created by the 
third pulse and detected during the third interval t^. The various pathways for the density 
matrix of a two level chromophore in Liouville space are shown in Fig 2. Correlations of 
the lineshapes during the flrst and the third interval provide information on environment 
dynamics during the intervening interval ^2- This supplements and greatly expands the 
information obtained from linear techniques. 



^analogous 2D NMR 
. Two-dimensional 
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43| . the picosecond 



2DCS can monitor dynamical processes at the femtosecond timesca' 
techniques are commonly used to study much slower (ms) processes 
infrared lineshapes have been used to probe is the structure of peptides 
hydrogen bonding dynamics by observing coherence transfer in molecular vibrations for 
phenol in benzene 4^ and for acetonitrile in methanol 45|. In the visible, 2DCS techniques 
have been used to study exciton transfer in photosynthetic antennae 46 1. 

Simulations of 2DCS signals usually employs either Markovian or Gaussian models for 
spectral fluctuations. The response functions for Markovian fluctuations may be obtained 
by the Green's function solution of the Stochastic Liouville equations 47, Q, [4^, which 



combine a Markovian master equation for jump dynamics with the Liouville equation for 
coherent evolution. The response functions of a multilevel chromophore linearly coupled 
to harmonic bath (Gaussian fluctuations) may be obtained by the second order cumulant 
expansion using the Wick theorem. All higher response functions may then be factorized 
into products of two point quantities. 
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In the present work we extend our earlier work [28|] to study signatures of aging in 2D 
lineshapes. In section II we build a general CTRW multistate jump model, and explain the 
condition of microscopic reversibility. The theory of 2D lineshapes is presented in section III. 
In section IV we discuss various parameter regimes of anomalous two state jump lineshapes. 
In section V we study aging effects in 2CDS spectroscopy and compare two approaches to 
describe aging: CTRW and time-dependent markovian master equation. The two models 
have same evolution of particle densities. The differences in 2DCS lineshape thus reflect 
the role of the underlying trajectory picture, i.e. unravelling of the master equation jl4|. in 
multipoint probes. 

II. TWO STATE CTRW JUMP MODEL; STATIONARITY, MICROSCOPIC RE- 
VERSIBILITY, AND AGING 



In this section we briefly review the anomalous relaxation model used in [2a|. The mul- 
tistate jump CTRW model is defined by a matrix ^(t) whose ij element is the waiting 
time probability density function (WTDF) for stochastic jumps from state j to state i. t 
is the time from the last jump where all memory is erased. The matrix is normalized as 



In the simplest two state jump (TSJ) model 48|, |49|, |50|, l5l|] bath has two states (a and 
h). We represent connection of density of renewals at various times in the a, h space by the 
matrix: 

\m ) 

The survival probability (j)i{t) (that no jump had occurred from state i for time t) defines the 
diagonal matrix of survival probabilities It is connected to the waiting time density 

function ipit) by (pjit) = [\I']-^- {t')dt' . The survival probability matrix thus connects 

the last renewal with final time 

The random walk is observed starting at time 0. The WTDF of the first jump ifj^t) may 
differ from 'ip{t) since it depends on how the system was prepared before t = 0. Similarly, 
the matrix represents the survival probability for the first jump. 



For a stationary process, the density of jumps to state i, rji is connected to the total 
density to be in state i, pi through 



pi = Piit) 



r],{t - t')<i)i{t')dt' 



poo poo 

Vi (j)i{t)dt = 'qi t^[^]ji{t)dt = rjiKi-i 
Jo Jo 



where we have used the fact that all densities p, rj are time independent for stationary process 
and Ki-i = t'^j[^]ji(t)dt is the mean waiting time in the i-th state. It then follows that 
all Ki-i must be finite. The rate for the j i jump is 

Jo 



We can now define the rate coefficients 



Aij = 



for i ^ j 



(3) 



The stationary density pj is thus obtained by the solution of the balance equation 

AijPj = 

j 

Using the same arguments, the WTDF for the first jump is 

Pj Jo "'I;* 



(4) 



The stationary condition (Eq.(jl]) is closely related to microscopic reversibility. CTRW is 
reversible if a trajectory ii, 12, ■ ■ ■ in with waiting times ^1,^2, ■ ■ ■ , (last time is survival) 
is equally probable its reverse in, ■ ■ ■ , ii with waiting times . . . , ^1. We thus require 

= [^'k-^iAU[%r.-2^r.-A^n^l) ■ ■ ■ (6)0n (6)P^„ (5) 

for all paths (sequences and waiting times). Eq. ([5]) can only be satisfied provided (i) the 
time profile of WTDF is independent of jump direction [\E']jj(t) = Tijipj{t) (5 2. l53l| . (ii) the 
rate coefficients for jump A^j = Tij/ni-j {i ^ j) must satisfy detailed balance |54 | 

^ijPj TjiPi 

1^1; j 
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and (iii) the probability of the first jump and the last survival are related through 
ip'i{t) = 4>i(t)/Ki-i which recovers Eq.(jl]). Eq. (jlj) thus expresses microscopic reversibility 
of a stationary ensemble: the survival probability coincides with the probability for the first 
jump backward. 

For the symmetric TSJ considered here (Eq. ([T])) we simply have 

and symmetric densities pa = Pb = ^f^. 

For < a < 1, Ki diverges, and it is impossible to construct a stationary ensemble. 



Asymptotically the jump rate decreases to l/ni which is in this case 3l|, l55|, ISGj. This 



scenario applies for arbitrary initial conditions. Many properties now depend on the initial 
observation time (aging). The normal diffusion constant for a Brownian particle moving 
on a lattice scales as ~ and its variance grows linearly with time (Einstein relation) 
(Ax^) ~ t/ni. When ni diverges, the particle loses its mobility at long times, and its 
variance (Ax^) growths grows sublinearly ~ t"' (anomalous diffusion). Another remarkable 
point is that the random walker survives at initial position for long times and ergodicity 
is broken. As a corollary, time averages obtained in single molecule measurements may be 
different from ensemble averages 34 1 . 

The simplest way to describe aging is by assuming that all random walks start by a jump 
made at some time to before the first laser pulse. The common choice ip'it) = V^(^) implies 
to = 0. The dependence on the initial observation time requires a to-dependent WTDF 
ip'{t;to). The consistent choice of ip'{t]tQ) will be discussed in section IV. 



III. SPECTRAL DIFFUSION IN 2DCS SIGNALS 

We consider a two- level chromophore with a ground l^f) and an excited state |e), transition 
frequency Qeg, and dipole moment peg subjected to three short laser pulses with an electric 
field E{t), and described by the Hamiltonian 

Hs = \e) [n,g + 6n,g{t)] (e| - E{t)p,g [\g){e\ + \e){g\] (7) 

5fleg{t) are stochastic frequency fluctuations caused by interaction with the environment 
and described by the CTRW dynamics. Observable quantities are obtained by averaging 
over all possible stochastic paths of 6Qeg{t) js?]. 
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We associate the frequency fluctuations with different bath states i, each inducing a 
transition frequency shift 6Qi. In TSJ the transition frequency assumes the value 6Qeg = 
(state a) and —Qq (state b). 

The response of our two level chromophore to three optical pulses is described by the 
third order response functions. The various contributions to the response function, known 
as Liouville space pathways (Fig 2), are labelled u [58|. During the intervals tj = tj — rj_i, 
between successive laser interactions the system's density matrix is in a given state |z/(^')) = 
|ee), \gg), \eg), or \ge) with corresponding frequencies qI^^ = 0, 0, Q,eg and —^eg respectively. 
The latter are modulated by the state of the bath. The Liouville operator describing the 
evolution in the bath state \eg); peg = L^gPeg is given by the following matrix in the a,b 
space 

Leg=\ " \, (8) 



iVlo 



where Lge = -Leg, and Lee = Lgg = 0. 



We next define the generating function p^ by the equation of motion. 

^ = -t5n,it)p. (9) 

with initial condition Pu{0) = 1. Here 6Qiy(t) = SQ^^\t) for t G (rj_i,rj). The third or- 
der response function for the i/'th pathway is then given by Rl (^3,^2,^1) = {Pu), where () 
implies averaging over the ensemble of bath paths. Coherent signals are generated only 
in specific phase-matching directions. Below we focus on the ki = — ki + k2 + ks and 
kii = ki — k2 + ks directions. In the rotating wave approximation these are represented 
by the four Liouville space pathways shown in Fig 2. 
The kj (photon echo) signal is 58 1 

yiih, h, h) = (^0 ;U^^e-^^-(*^-*^) [Ruih, h, h) + RUh. t2, h)] , (10) 
and the kii signal is 

yiiits, t2, h) = (^0 /x^^e-'^-(*i+*3) [Riih, t2, ti) + Riiiits, t2, h)] (11) 

For stochastic models such as considered here the bath evolution and equilibrium state are 
independent on the state of the system pee or pgg so that Ri = Rm, Ru = Ri^. 
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The third order correlation function for the i/'th pathway may be obtained by solving 
Eq. © 



Rf\h,t2,h)=e{h)9{t2)9{ti)le^Y> 



T3 







exp 






L TO 



:i2) 



where the upper sign represents Ri = Rm and the lower Ra = Ri^. 

The 2D signals are defined by frequency- frequency (ci;3,a;i) correlation plots for a fixed 

t2. 




(13) 



We shall also display the following combination, which shows simpler lineshapes with 
purely absorptive peaks [59|, [60 1 . 



5'a(c^3, t2, UJl) = Sl{uJ3, t2, -UJl) + Sll{uJs, t2, UJl) 



(15) 



The response is represented in a, b space by a matrix whose jl element accounts for 

(3) 

the contribution to Rl from paths with an initial bath state / and final state j. 



{t^,t„t,) = J2 [G% ih, hM) [pAi {t = 0) 
jl 



(16) 



For Markovian relaxation [\I']jj(t) = Tije~*f'^^'^ / Ki-j each G'^ may be factorized into a 
product of three Green's functions representing the time evolution during the ti ,t2 and ts 
intervals whereby the density matrix is in the and u^^^ states. 



^^^(^3,^2,^1) = G''''\ts)G''''\t2)G''''\h) 



(17) 



The Green's functions can be calculated by solving the stochastic Liouville equations (SLE) 
47|. 

where A is the matrix of jump rate coefficients (Eq. (JSj)). The SLE has recently been applied 
to describe vibrational 2D signals for frequency fluctuations modulated by hydrogen bonding 
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of phenol in benzene 6l| , conformation changes of peptides 6J] and infrared hneshapes of 



water 



63|. 



The simulation of systems with long memory is much more complex. Various types 
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20| for 



of reduced equations of motion for the CTRW dynamics have been developed 
calculating the two-point correlation functions. These, however, may not be extended to 
multipoint quantities required for the description of 2DCS |6J], since the factorization, Eq. 
ffTTll. does not hold for nonmarkovian relaxation. 



We have recently 6J] developed an algorithm for solving this model. This is based 
on the successive recurrent construction of a hierarchy of Green's functions. It relies on 
the renewal property computing the CTRW [57]. Below we present an alternative, more 
intuitive, derivation which is reminiscent of the Green's function method. 

We need to maintain a bookkeeping of whether or not there was a jump during each of 
the three time intervals ti, t2, t^. For each of the three intervals we must distinguish between 
two possibilities; either there was no jump or there was a least one jump. is thus given 
by a sum of 2'^ = 8 terms each representing one type of path in bath space. 



G^ih^hyti) — G'^jts, t2, ti) 



(18) 



m=l 



These terms are depicted in Fig 3, where the presence of any (> 1) jump in a given time 
interval is represented by the trajectory touching the time axis. 



G'^ are conveniently recast in Laplace space. We define (our notation is similar to |65| ) 

^{s-L) = e-''^{t) exp (^Lt^ dt 
This implies for our TSJ model 



eg J 



i){s-ino) 

V'(s + ifio) 



(19) 



where ip{^) = Jq° i^it)^ '^^dt is the Laplace transform of ip. ^{s — L) for the survival function 
is defined similarly 

~)(s + i^o) 

4){s - iQo) 

G^ is expressed as a matrix product of the propagators through the intervals with any 
jump in the particular interval ( if the trajectory touches the axis in Fig 3) , with additional 



$(S - L^g) 



(20) 
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factors for segments connecting different intervals. These ensure that the bath state does not 
change between the last jump in the earlier interval and the first jump in the later interval. 
Both factors will be described below. 

We first calculate the evolution for a fixed state of bath where no jump occurs over several 
time intervals. Let us assume that the state is fixed for time t'^ in the m-th interval, till time 
t'l in some subsequent 1-th interval and during all the intermediate intervals ti, I > i > m. 
The probability of this evolution is either or $' depending on the path. The 

propagator connecting the state immediately after and after t'^ is given by 

t(vi/, t;, t,_,, . . . , 4) = ^(t; + 4+5^ exp L(')t; + l(-)4 + ^^^^^ (21) 

i=m+l \ i=m+l / 

where L*^*^ = ±Leg, 0, depending on the state of the density matrix in the i-th interval. This 
contribution may appear in several ways. Either for the evolution between the last jump in 
the m-th interval and the successive jump, first in the 1-th interval, or for the very first jump 
when the interval does not exist and \1' It also appears for the survival from the 

very last jump when disappear and \1/ ^ $. Finally, when no jump occurs, then \1/ 
t'l is absent and t'^ = ti . 

The second ingredient in our calculation is the propagator through the k'th interval 
described by the integral equation 

S(r) = [ ^(r - r') exp \-iL^''\t - r')] t{T')dT' (22) 

with S(0) = i. The matrix S connects the arrival densities at two times within the same 
interval tj. By solving Eq. (!22l) in Laplace space, we obtain the following propagator through 
the k-th interval 

t{s,)= (23) 

This contribution appears provided some jump had occurred in the k-th interval, (i.e. the 
trajectory touches the axis in the k'th interval in Fig 3.) Eq. (1231) can be interpreted as a 
summation of a geometric series for paths with 1,2, . . .jumps in Laplace space, where time 
convolutions become simple multiplications. 

All of these factors should be convoluted in time to generate the trajectory. For instance, 
the domain of integration for the first contribution Gi is shown in Fig 4: 

GUt3,h,t,)= di, di, diJ d^J di2 di^ 

Jo Jo Jo Jo Jo Jo 
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X t($, ^6)2(^3 - ^6 - e5)t(vl>, ^5, ^4)t{t2 - ^4 - ^s) t(vl>, ^3, 6)S(tl " 6 " 6)t(vl/', ^l) (24) 

This results in a simple product in Laplace space 

C(^3, ^2, Si) = t ($, S3)S(s3)'f (Vj/, S3, S2)S(^52)'f ^2, Si)S(si)t (v]/', Si) (25) 

We have already calculated Laplace domain S (Eq. (!23|) ). T can be easily transformed 



as well, leading to equivalent results to those reported Appendix C of Ref.[6J]. Eq. fl2Sl) is 
finally expanded in terms of the matrices $, ^ and the complete expressions agrees with 
Appendix D of Ref.[64], where was obtained in a different way. 

Since the response functions (Eq. f|T2|) ) are causal, the 2D lineshapes (Eq. [T4|) may be 
obtained by analytical continuation of si,S3 (the Laplace variable conjugate to ti and ^3). 
The ^2 variable is obtained by reverse Laplace transform using Bromwich integral 

Sl{uJ3, t2, -(^1) = ^-^"^ J dS2e^^*^Rii (S3 = -Z(CJ3 - Vteg), ^2, Si = - Vteg)) 

^4 /•too 

Sll{uJz,t2,UJi) = -^Im / dS2e^^^^ Ri{s3 = -i {UJ^ - Qeg), S2, Si = -i{uJi -^eg)) 

J -ioo 

For t2 = these integrals may be calculated analytically. The resulting two-interval 
functions may be alternatively obtained by directly building the two interval (t3,ti) response 
function. 

IV. LINESHAPES FOR STATIONARY ANOMALOUS RANDOM WALKS 

Microscopic reversibility in stationary ensembles implies that ^i{t3,t2,ti) = 
—S^j'{ti,t2,ts) which in the frequency domain gives 

Si{uJs, t2, -uji) = Si{uJi,t2, -ujs) (26) 

Similarly =5^/7(^3,^2,^1) = -^7/(^1,^2,^3) which implies 

Sii{uj3, t2, uji) = Sn{ui, t2, UJ3). (27) 

Combining Eqs fl26p and fl27p with Eq. fllSp we obtain the following symmetry of the 
lineshape 

Sa{uJ3, t2, UJi) = Sa{uJi, t2, UJ3) (28) 
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Thus 5"/, Sji, and 5"^ are symmetric to the interchange of uji and 003. 

When during the ^3 interval the bath has lost its memory of its state during ti (e.g. 
normal relaxation with t2 00), the response functions may be factorized as 

yj{t,, h, h) = 2{i/h)K{t,)K*{h) (29) 

and 

^//(ts, t2, ti) = 2ii/h)Kits)Kih) (30) 

Here 

K{t) ^ (z/;i)/.^^e-^-*(exp[-z / 5fie,(r)rfr]) 

is the linear response function for stationary ensembles. Its Fourier transform gives the 
absorption lineshape 

POO 

WAiuj) = Im K{t)exp[iujt]dt (31) 
Jo 

(The absorption of a nonstationary ensemble is not proportional to the Fourier transform 
of the linear response function [50].) 

Using Eos. (129|) . (!30l) . and (!3T!) . Sa then reduces to the product of the linear absorption 



lineshapes 



66| 



hSAitOs, t2 ^ 00, CUi) = AWA{uJi)WAiuJs). (32) 



Algebraic memory decays will result in a slow convergence to this asymptotic lineshape. 
In addition, as will be shown below, the spectra diverge at certain frequencies where the 
factorization (Eq. (l32l) ) does not hold. 



We shall consider a specific model of anomalous relaxation with the WTDF 30|, IfiJ]: 



V'(s) = TTT—, ^ — tt; 1 < a < 2; (33) 

Ki is the mean of '?/'(t), while ka controls the long time algebraic tails il)w{t) ~ 
Note that Eq. ([6]) may be conveniently represented in Laplace space 

SKi 

We first consider the linear response obtained from the one-interval Green's function 
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The kernel may be calculated by 



Q(t) = t($',t)+ fd^2 f ^'rf6t($,6)s(ti-6-6)t(^',ei) 



Transforming into the Laplace space domain yields 



50 



64 



65] 



-1 



g(s) = $ (s-L) + <l>(s-L) l-^(s-L) m{s-L) 
Combining Eqs. fl3T|) . fl33|) . and fl34l) we finally get 



(34) 



X i?e 



1 



(35) 

In all plots we use dimensionless frequency units {ijjj — VLeg)/^o by setting VL^g = 0, fio = 1- 
In Fig 5 we display the absorption spectrum in the slow (/ti^^o > 1? top) and the fast 
(ki^^o < 1) bottom) fluctuation limits. The lineshape has two peaks at = ±1 and in the 
fast fluctuation limit we obtain a flnite central peak 50|, l6^ • The fraction of particles that 
remained at the initial position is signiflcant (not exponentially small) at all times. This 
results in the divergence of peaks at = ±1 



W{uj) ^ C\Auj 



ia-2 



(36) 



with C = cos [7r(l — q;/2)]k^ /2, and where t 



and Au = u — Qeg + for uj = —1 peak 



5G 



le detuning is Aa; = u — Q^g — ioi u = 1 



64| . 



The parameter a controls the peak singularity. For a — 2 the divergence is cured and we 
approach the Markovian lineshape. For fast fluctuations ^0^1 << 1 the central peak grows, 
as Ki becomes shorter. This is reminiscent of the motional narrowing for the Markovian 
case. However, the two divergent peaks still retain an anomalous lineshape. 

In Fig. 6A we display the Sj^u^, —Ui), Sji{u!3,u!i) and S'^(co'3, u^l) signals for slow fluc- 



tuations QqKi » 1 and ^2 = 0. Similar to the Markovian case [61i, all panels show two 
diagonal-peaks at (^73,^^1) = (1, 1) and (—1, —1). However the peaks are nonlorentzian and 
divergent. Sj and Su diverge along the ui = ±1, = ±1 lines, but much of this divergence 
is cancelled in 5*^ which only diverges at peaks (1,1), and (-1,-1). 

We next examine the analytic structure of these divergencies for the ui = 1 and 00^ = ! 
lines. The slowest decay is connected with the survival function for the flrst jump 0'(t) ~ 
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t^~°'. Gf, is thus the most rapidly divergent term. The analysis of peak divergencies thus 
reduces to the Gl contribution. We denote Atus = — Q^g — and Aui = ui — Qeg — 
and find 

Sl8{(^3, t2 = 0, -UJi) = -T^Im — 

Snsiu^s, h = 0, 0.0 = ^Im _ (37) 

The lineshapes (Eqs.(39)) diverge along the lines Au^ = and Aui = 0. The divergent peak 
structure is summarized in Table I. The left column corresponds to situation when Auj^ is 
held fixed at a small but nonzero value and Aui approaches the singular point 0. Thus we 
consider Auji << Au^ and (f)'{i Aui) » cj)' {—iAuo^). With the asymptotic expansion 

~ (38) 

we get the asymptotic form of divergent Siiuj^, — cui), and Snioj^^oJi) shown in the Table I. 
In the right column we similarly approach the singular line at AcJa = 0. We have verified 
these analytic asymptotic results numerically (not shown, it also qualitatively agrees with 
Fig 6A. 

Si and Sn have opposite signs, and their combination Sa is finite due to interference. The 
divergencies are only seen at the (1,1) and (-1,-1) peaks, and not along the entire Ui = ±1 
and LJ^ = ±1 lines, since the Sj, Sjj divergencies cancel. 5*^ is finite, but nondifferentiable 
along these lines. 

We next examine more closely the variation along the AuJi = axis. 

SA8{AiU3, t2 = 0, AuJi = 0) = 

The asymptotic expansion Eq. fl38l) yields the analytical peak structure at Auj^ ~ 0. 

Sa{Auj3, t2 = 0, Atoi = 0) ^ BAco^-^ (39) 

5 = ^/t^-isin [7r(2-a)/2] 

The analytic structure of the (-1,-1) peaks is the same. This follows from the assumed 
[^]ab = [^]ba symmetry of TSJ model, which implies 5',^(ti;3 + ^eg, "^i =F ^eg) = 'S'i^(— ^^3 + 
^eg, h, — t^i =F ^eg)] Upper sigu applies for z/ = J lower for z/ = II, A. Based on Fig 6a, the 
peaks are more localized with steeper contours for smaller a. In all cases we see a dip at 
(0,0). 
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The two peaks induced by Gg are universal and survive even for the case of fast fluc- 
tuations QqKi « 1, as shown at Fig 6B. Rapid changes during ti and induce a new 
peak at the average frequency (0,0), (motional narrowing). The Sa_ (0,0) peak is Lorentzian: 
The star-like contours, best seen for a = 1.2 correspond to a product two Lorentzians along 
uji and u^. The Sj and Sn lineshapes are similar. Both may be described by a statisti- 
cal mixture of rapidly fluctuating particles responsible for the central peak, with the static 
phase responsible for the divergent peaks at the fundamental frequencies. Surprisingly, this 
picture is most pronounced for small a = 1.2 where all peaks are well-separated. Increasing 
a broadens the (-1,-1), and (1,1) peaks, making them interfere with the central peak, and 
the lorentzian shape becomes less pronounced as a — > 2. 

The variation of Sa with t2 in the slow fluctuation limit is displayed in Fig 7. For t2 
longer than the mean waiting time ki fractions of trajectories have different frequencies in 
the ti and ts intervals , as described by the Gq contribution resulting in new cross peaks at 
(-1,1); (1,-1). Since we are in the slow fluctuation limit the peaks are still well resolved. Both 
diagonal and cross peak contours are elongated along the ui^^ = ±1 directions. Nevertheless 
the decay of the Gq contribution t^g^ (compared to the diverging t^g^ decay of Gs which is 
relevant for diagonal peaks) is integrable and thus the cross peaks do not diverge. Another 
notable point is the breakdown of Eq.( l32i) at cui^a = ±1; Memory loss is not complete since 
the algebraic functions do not factorize. At other frecjuencies the lineshapes approach this 
limiting lineshapes (Eq. (!32l) ) algebraically as " [28]. These simulations illustrate the 
capacity of 2DCS to probe anomalous relaxation during the t2 interval. 



V. NONSTATIONARY ENSEMBLES; AGING OF 2D LINESHAPES 



In our earlier work [281] we considered nonstationary ensembles with < a < 1 by 
assuming that the random walk is started by a jump at the time origin, coinciding with the 
first laser pulse, so that response may be calculated hj ip' = 4'- The lack of microscopic 
reversibility is reflected in violations of the symmetry relations Eq. ( 128|1 . The higher mobility 
during the (earher) ti interval compared to resulted in broader peaks along the Ui axis 
compared to u^. 

Here we explore signatures of aging. We consider random walks, which start by a jump 
made at some time to before the first laser pulse and examine how the nonlinear lineshapes 
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vary with to- The response function then depends to even for to ^ oo. This is known as 
aging. All aging effects are fully described by calculating the WTDF \E''(t;to) for the first 
jump which is now to dependent which must be consistent with the CTRW dynamics during 
the to period. 

\E''(t; to) can be calculated along the lines of Eq. (121 p by omitting the coherence evolution 
= during to, 

^'(t,to) = ^(t + to)+ Trf^ r + 6)2(^0 -6-6)^(^1) 

Jo Jo 

In Laplace space we find for our TSJ model 

^'{s-so) = ^M^^^ (40) 

[1 - ^'(■5o)](-5 - So) 

The 2D lineshapes may thus be calculated using the algorithm presented in Section III. The 
to dependence is obtained by numerically inverting these Laplace domain formulas. 

The long to limit may be obtained by setting so — >■ 0. For CTRW with finite Hi the 
denominator in Eq. ( I40l) is 



1 - V(so) ^ -so- 



— f^iSo 

s=0 



ds 

This reproduces the WTDF of the first jump for a stationary random walk tp'{t] t^) = (l){t)/ni. 



The lack of stationarity has some important consequences. As pointed in 501] frequency 
domain absorption measurement is no longer given by the Fourier transformed response 
function. Thus the absorption of an aging ensemble can not be calculated using Eq. (1341) . 
Fortunately, 2DCS works in the time domain, and the measurement directly probes the 
response function. Thus the problems discussed in [50] do not apply for impulsive time- 
domain techniques such as 2DCS lineshape. 

A more subtle point is that due to the lack of equilibration, the averaging over consecutive 
pulse sequences may depend on the experimental data acquisition repetition rate. Proper 
definition of the response function requires a careful preparation ip'{t) before each pulse 
sequence. 

We have calculated the response functions the variation of the lineshape with the prepa- 
ration time to for the following model 

^^(s) = ^^^, «G(0,1) (41) 
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This corresponds to a WTDF with algebraic tails iplt) ~ (/t/t)^"'"". 

We took a = 0.98, which is close to the Markovian case (Eq. fHTj) for a = 1) in the 
fast fluctuation limit kQq << 1. This choice is motivated by the simpler interpretation of 
the lineshapes; we expect it to be closer to the Markovian case than the rather complex 
to = shapes presented in [28|]. The effect of to could thus be better isolated. In addition, 
aging effects appear at arbitrarily long timescales (for suitable choice of parameters). This 
overcomes the difficulty with strongly anomalous ensembles, whose lineshapes cannot be 
obtained by repeated measurements on the same sample, whose response function is changed 
between two pulse sequences. 

The top left panel of Fig. 8 (to = 0) shows fast-fluctuation Markovian contours and only 
tiny peaks at (1,1), and (-1,-1). No signatures of time irreversibility are seen since Eq. f l28p 
is nearly satisfied. We next increase the aging time to as we move from the top left panel to 
bottom right panel. The (1,1) and (-1,-1) peaks appear and grow, while the central peaks 
slowly get weaker. This reflects decrease of the jump rate with time. Some small deviation 
from the symmetry relation Eq. (|28l) can be noticed. The process is nearly reversible on 
the ^0 ^ timescale which dominates the lineshapes. A remarkable point is that the central 
(motional narrowing) peak coexists with these static limit peaks. The anomalous process is 
better viewed as a mixture of static and fluctuating particles, rather than a homogeneous 
rate. 

This clearly distinguishes our algorithm from calculations based on time-dependent rate 
master equations, which do not allow to properly describe memory effects in multipoint 
probes. To support this statement we have constructed Markovian process subjected to the 
same master equation, i.e. we require correct prediction of total densities and subsequently 
apply them to calculating response or multipoint correlation function based on Markovian 
schemes. The trajectory picture of both approaches is different [g^]. 

Consider a Markovian master equation whereby densities evolve in the same way as the 
aging random walk for arbitrary initial densities, i.e. it has the same Green's function G{t) 



The master equation is constructed by differentiating Eq. (H2l) with respect to time 



Pit) = G(t)p(0) 



(42) 



dp{t) 
dt 



A(t)p(t); 



dGit) 



dt 



G-\t). 



(43) 
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The transition matrix A of time-convolutionless master equation is thus uniquely defined. 
The Green's function Eq.f l42p is the solution of the mater equation 

G{t) = expT [ A{t')dt' (44) 
Jo 

We consider a symmetric two state dynamics parametrized by a single function A 



Ait) 



Eq. (l44l) can be solved after a simple algebra. This gives 

G'ii(t) - Gio(t) = exp [-2 [ A{t')dt'] (45) 

Jo 

Inverting Eq. fl45l) . the rates can be calculated once the Green's function is known 



. ^ -||GuW-G.o(i)| 

2|G„(i)-G,„(i)l ' ' 

We next adjust the Green's function to agree with those of our aging random walk. In 
Laplace space it reads 

r (,] r ( - - A^Jti^ 

^U[S - ^io[S — — — -r^ — -J— — -r-TT 

l + tp{s) s[l + 'ip{s)\ 

For the model Eq. (14T!) . 

Gu(s)-GMs) ^""'^^ 



s[l + (ks)"] 

which may be also calculated directly in time domain as series 



Gn{t) - G,o{t) = ^(-1)" X (47) 
^ r(na + 1) 

with the gamma function T{y) = x^^^e'^dx. The master equation is thus defined 
by combining Eqs. (H3|) . P6|) . and (147|) . The rates decay asymptotically (t oo) as 
A{t) ^ a/{2t). Exponential WTDF's {a = 1) correspond to constant rate A = 

Fig 9. shows the time dependent rate of the master equation for various a. Aging 
effects (decreasing mobility with time) are refiected in the decreasing rates. Increasing a the 
decay is slower when approaching the markovian limit (a = 1) and the rates change slowly 
for long periods. This regime is particularly interesting because it may provide sufficient 
time to measure the rate constant by e.g. lineshape experiments and give clear meaning 
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to our arguments. (Diverging rates at very small times t << k are integrable and thus 
insignificant.) 

We shall compare two types of stochastic processes subjected to the same master equa- 
tion, but with different unravelling into trajectories Aging lineshapes for the CTRW 
model were already presented at Fig. 8. The second model is defined by Markovian pre- 
scription: The probability of jumps are independent of the past trajectory. The stochastic 
Liouville equations and Green's function technique may then be used to calculate the non- 
linear response. 



RAh, t2, ti; to) = 0{t3)9{t2)9{ti)( expT 



expj, 



n 



L-J TO 



We are interested in the peak pattern, which is influenced by fluctuations on the 
timescale. We consider a parameter regime where the rate does not change significantly on 
this timescale, and thus the peak pattern may be analyzed by a simple approximation of 
rates independent of ti, and and analyze aging of ^2 = lineshapes 



R^ih, 0, h; to) = ^(t3)^(t2)^^(ti) (exp [(^(to) + L^^)) t 



exp 



A(to)±L«)ti 



We then obtain 



- ujscoi =F ^0 ^ i2A{uJs + coi] 

-He- 



(48) 



[ul - 1^0 + ^2cJiA] [cj| - 1^0 + «2cj3A] 
where the upper (lower) sign is for {3 = I {(3 = II) and where A = A(to). 
The aging Markovian 2D absorptive lineshapes are presented at Fig 10. The central peak 
is gradually broadened with increasing time (and decreasing rates) and splits into two peaks 
centered along diagonal at fundamental frequency. These peaks get narrower for long t2. 

The significance of the different trajectory picture can be seen by comparing the two 
lineshapes at Fig 8 and 10. We notice that the crossover to static lineshapes is somewhat 
faster at Fig 10. This is, however, less obvious feature, since it depends on chosen particular 
parametrization. The more significant feature, which distinguishes the two models is that 
the static peaks at fundamental frequencies and the fast motional narrowing central peak 
never coexist at Fig 10 in contrast to Fig 8. 

This may be explained as the direct signature of memory. The CTRW model shows 
two populations static and fast fluctuating, i.e. particles are differentiated based on their 
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histories. In contrast all particles in the Markovian model have homogenous probabilities 
for the next jump. This lack of memory is reflected in the unique peak pattern with no 
simultaneous static and fast fluctuating signatures in the spectrum. 

These two models are nonequivalent since they assign different trajectory picture to the 
same density matrix, as is clearly seen from the higher order correlation functions and 
response. The coexistence of both static and fast fluctuations in spectra clearly reflects the 
additional information, beyond the two point correlation functions. 

The unravelling of master equations into trajectories is an important issue. Two dimen- 
sional lineshapes which are sensitive to the trajectories should provide a direct test for the 
unravelling schemes Single molecule spectroscopy looks at the trajectories one at a 

time. Multidimensional spectroscopy looks at the entire ensemble but unravels it by the 
manipulation of coherence. 

In summary, our simulations demonstrate that two-dimensional correlation plots of signals 
obtained from the response of the system to sequences of multiple laser pulses carry specific 
and direct signatures of complex dynamics. Such techniques are currently feasible in many 
spectral regimes, NMR, EPR, the infrared (vibrations, phonons) and in the visible (electronic 
excitations). 
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Table captions 

Table 1 Sj, Sn lineshapes shows divergent growth along the cu^ — D,eg—^o and cui — ileg—fto 
lines. Table I shows their asymptotic form. 
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Figure captions 



Fig 1 Pulse configuration and time variables for a four wave mixing experiment . 

Fig 2 Feynman diagrams for the third order response of a two level chromophore with 
wavevector ki = — ki + k2 + ks and kn = ki — k2 + ks. 

Fig 3 The 8 contributions to the Green's function (Eq. f|T8|) ) of the third order response. 
Contributions represent paths with (when the line touch the axis) or without (when 
the line does not touch the axis) some jump during each of the three time intervals 

^1,^2,^3- 

Fig 4 Integration time variables in Eq. ( 12^ . 

Fig 5 (Color online) Linear absorption for slow kiQq = 2 (top panel) and fast kiQq = 0.2 
fluctuations and different a as indicated. ^IqKa = 0.5. 

Fig 6A (Color Online) The Si{u!3,0, —ui) (top), Sjj{u!3,0,uji) (middle), and Sa{uj3,0,uji) 
(bottom) signals (Eq.([T5])) for the WTDF (Eq. ([33])) at = 0, for slow fluctuations 
QoKi = 2.0, and k^/ki = 0.25, a = 1.2 (left), 1.5 (middle), 1.8 (right). 

Fig 6B (Color Online) The same as Fig 6a but for fast fluctuations QqKi = 0.2, and 
ka^o = 0.5, a = 1.2 (left), 1.5 (middle), 1.8 (right). 

Fig 7 (Color Online) The Sa signal (Eq.([l5])) for the WTDF (Eq. ([33])) for (left to right) 
a = 1.2,1.5,1.8, and ka/ = 0.25, QqKi = 2.0, ^2 = i^i (top), t2 = 2Ki(middle), 
^2 = 10/ti (bottom). 

Fig 8 (Color Online) Aging effects in 2D lineshapes. The Sa signal (Eq. (fT5]) ) for the 
nonstationary random walk model (Eq. fl4Tl) ) ^2 = 0, kQq = 0.2, a = 0.98 for various 
initial time (from left top, to right bottom) to = Ok, 10k, IO^k, IO^k, IO^k, IO^k. 

Fig 9 (Color online) Time-dependent rates of aging random walk Eq. fll6]) for a = 0.3 
(sohd), 0.5 (dashed), 0.7 (short-dashed), and 0.98 (dotted line). 
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10 (Color online) Aging 2D Markovian lineshape (Eq. (HHIl ) for various initial time 
to = 10~^, 1, 5, 10, 20, and 100 . Master equation for probability densities correspond 
to the random walk of Fig. 8 
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fixed 




Aa;i 


varied 




Aa;3 


Siius, -ui) ~ ^ sin [7r(2 - a)/2] x 




-sgn(Ac.3)IX"^ 


Sii{uJ3, uJi) ~ ^ sin [7r(2 - q;)/2] x 




sgn(Ac.3)^ 



Table I 
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Fig 5 
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